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Abstract 

We study a two-dimensional granular system where external driving force 
is applied to each particle in the system in such a way that the system is driven 
into a steady state by balancing the energy input and the dissipation due to in- 
elastic collision between particles. The velocities of the particles in the steady 
state satisfy the Maxwellian distribution. We measure the density-density 
correlation and the velocity-velocity correlation functions in the steady state 
and find that they are of power-law scaling forms. The locations of collision 
events are observed to be time-correlated and such a correlation is described 
by another power-law form. We also find that the dissipated energy obeys a 
power-law distribution. These results indicate that the system evolves into a 
critical state where there are neither characteristic spatial nor temporal scales 
in the correlation functions. A test particle exhibits an anomalous diffusion 
which is apparently similar to the Richardson law in a three-dimensional 
turbulent flow. 
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I. INTRODUCTION 



Granular media have distinct behaviors from those of usual solids, fluids and gases . 
There are two particularly important features that contribute to the unique properties of 
granular materials: ordinary thermal fluctuations play no role because of the mesoscopic 
size of grains, and the interactions between grains are dissipative 0. Therefore, in order 
to maintain the dynamics of granular media in the long run, external driving forces are 
inevitable. Both experiments and computer simulations show that the dynamical responses 
of granular media to external forces exhibit a wide variety of interesting phenomena |T] fj . 
Traditionally, external energy flows into the granular system from the boundaries, either by 
shear or by other means [0-fJ. 

The behaviors of inelastic and elastic systems under the same situation are quite different. 
Elastic systems, if contacted with a heat bath where energy flows into the system from the 
boundaries, will attain a homogeneous state with a uniform temperature (average kinetic 
energy) field. However, the local kinetic energy of granular systems exhibits spatial gradients 
due to inelastic collision if energy is input from the boundaries . A freely evolving granular 
medium may have uniform temperature, but that temperature is time-dependent [|irj.|Tl] . 



A uniform and time-independent temperature is essentially required by the granular 
thermodynamics theories |T2}^H| or by any effort to make a comprehensive comparison be- 



tween granular materials and fluids at equilibrium. Recently there are investigations []lq-|l7 



showing similarity between the steady states of granular materials and usual fluids at equi- 
librium. 

Thus in order to study the properties of granular media free from the nonuniformity of 
temperature and anisotropy due to gravity, one needs to introduce an idealized system. In 
this paper, we carry out computer simulations in two dimensions of a model system satisfing 
the above conditions. The driving force is applied to each particle so that the external 
energy flows into the system uniformly. 

Our main concern is the correlations in the steady state of the granular system where 
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energy dissipation due to inelastic collision is balanced with the energy input. The spatial 
distribution of collision events, particle diffusion, density-density correlation and velocity- 
velocity correlation will be investigated by simulations. The results show that there are 
no characteristic spatial and temporal scales in the correlation functions, no characteristic 
energy scale in dissipation, and an anomalous diffusion of a test particle. 

We admit that the way of energy input uniformly to each particle in our model is difficult 
to realize in experiments. As we mentioned above, however, this kind of idealization is 
necessary to explore the properties inherent to systems with inelastic interations and to 
clarify the essential difference between granular materials and usual fluids. 

This paper is organized as follows. In Section II we present the model. Numerical results 
from the model are given in Section III and we close the paper with Section IV which is 
devoted to discussion. 



We consider N hard disks (particles) of diameter a interacting via collisions in a two- 
dimensional square cell of length L with periodic boundary conditions in the x- and y- 
directions. During a collision, the relative velocity of two particles in the normal direction is 
reduced by the restitution coefficient R while in the tangent direction the collision is elastic: 



where primes indicate velocities after the collision and n is a unit vector pointing from the 
center of ith particle to that of jth particle. The loss of energy due to the collision is given 



II. THE MODEL 



v; = v, - i(l + R)[n -(Vi-v^fi, 



(1) 




v^ + Kl + i^Mvi-v^n 



by 



^- = im(l-i2 2 )[n.(v i -v i )] 2 



(2) 



where m is the mass of the particles. 
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On the other hand, energy is input to the system after each time period of T in the 
following way: the velocities are perturbed instantaneously by a random amount: 

< = v, + & (3) 

where primes refer to the velocities after energy input. The components of £ x and £ y , 
are random numbers distributed uniformly between —SV and SV. The energy input rate is 
therefore proportional to (SV) 2 /T. The way of adding energy to the system in Eq. @ can 
remove any macroscopic flow of the system since there is no preferable direction chosen. 

One-dimensional version of the present model was studied by Williams and MacKintosh 
18] with continuous energy input. Here the energy is added to the system instantaneously 



after each time period of T. Therefore, we are able to use the event-driven algorithm |T9|,[7[| 
to perform the simulations. An event is defined as a change of velocity either by collision or 
by energy input. 

The inelastic collapse singularity was observed in simulations of inelastic system |20| . 



This is caused by the appearance of infinite number of collisions between a group of particles 
which are spontaneously arranged on a straight line [ZIJ. We avoid this singularity by a 



slight modified collision rule as introduced by Deltour and Barrat [ ID | . At each collision, the 
relative velocity of the two particles is first computed (Eq. ([[])) and then rotated randomly 
by less than 5 degrees JT0| . 



The parameters in our systems are the restitution coefficient R, number of particles N, 
system size L and average density (area coverage ratio) p. The diameter a of particles is 
related to them by p = . We will keep L = 1 throughout the simulations, and show 
the results for SV = 0.001 and T = 0.002 unless stated otherwise. The heating period of 
T = 0.002 corresponds to adding energy to the system about every five collision events for 
# = 0.1. 



III. NUMERICAL RESULTS 



4 



A. Approach to steady state 



Starting from any initial configuration the system with dissipation (R < 1) reaches a 
steady state in the center of mass frame. Since the energy input (Eq.(|3])) does not guarantee 
conservation of the total momentum, we subtracted the average velocity from the velocity 
of each particle in order to remove the motion of the center of mass of the whole system. 
Thus the velocity is defined in the frame moving with the center of mass as it was done in 



Ref. [18 



Fig. 1 shows how the average energy per particle relaxes with time towards a steady 
state for a system with N = 1024, p = 0.16, and R = 0.8. It is found that the relaxation 
is exponential as E(t) = E(oo) + (E(0) — i?(oo))e~*/' r . We find that the relaxation time 
r (= 2.63 in Fig. 1) is independent of the initial state and is controlled by the model 
parameters. 

Theoretically, we can derive the energy relaxation equation as follows. The system 
increases its energy as a result of external driving (Eq.(|3|)) with an amount of 

-> N i N i N 

A£ + = \ £(v, + £J 2 -\Y.< = \ + 2v, ■ Q (4) 

i=l i=l i=l 

while it decreases its energy at a collision according to Eq.(^J). The energy change rate is 

AE 



aAE + - bAE- (5) 



At 

where a and b are constants, depending on the frequency of energy input, number of particles 
and collision rate. 

After averaging over a time interval much longer than the collision time (which will be 
defined later), Eq. (^) becomes 

£ = <t-VE (6) 

where a' =< |a£ 2 >= ^a(5V) 2 , b' = \m{l — R?)b and we have made use of < Vj • £j >= 0. 
A steady state with E(oo) = a'/b' exists for Eq. (|5[). The time-dependent relaxation is 
given by 



E(t) = E(oo) + Ae- Vt 
where A is an integral constant depending on the initial condition. 



(7) 



B. Collision rate 

In the steady state the collision rate is constant, independent of time. This is in con- 
strast to the time dependent collision rate in the freely evolving granular medium where the 



total number of collisions increases as ln(l + t/t e ) for homogeneous cooling [jTQ] , p~Tf| . In our 
simulations, we count the total number (C) of collisions starting from a configuration in 
the steady state. Fig. 2(a) shows that the total collision number is linear with time for all 
systems with R < 1. The slopes of the lines are the collision rates for different values of R. 
One would expect that the collision rate should be an increasing function with increasing 
R as the velocity (therefore granular temperature defined below) is larger and the collision 
is more likely to happen in the larger R system. However, we find from Fig. 2(a) that the 
collision rate increases as the degree of dissipation increases. Fig. 2(b) is the dependence 
of the collision rate 7 as a function of degree of dissipation e = 1 — R 2 . We find that the 
numerical data can be fitted using the following form 

where the exponent c is found to be about 0.63 for the data in Fig. 2. It is evident that the 
collision rate will diverge if e is equal to 1. 

The increase of collision rate with decreasing R may be linked with the clustering mech- 
anism found in Refs. p,2l|. However, by direct visualization the density nonuniformity is 



not as clear as in the situation of the freely evolving granular medium @,|2T]]. Fig. 3(a) is 
a typical configuration taken from the simulation. To find out any correlation in the col- 
lision events, we therefore resort to the following method. We record the positions where 
dissipation takes place by, e.g., taking the middle point's coordinates between each colliding 
pair. At each dissipation place, we draw a data point and by doing so for some time interval 
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we obtain Fig. 3(b). We see clearly that there is clustering in the dissipation places. Note 
that the clustering is not due to limited data points (or short time observation). In fact, 
there are 20,480 points in Fig. 3(b) with overlap. We have confirmed this observation by 
longer runs. It is also verified that such a clustering phenomenon does not exist for elastic 
system without dissipation (R — 1). We observe that the clustering patterns change their 
positions and forms as time develops. The clustering will be quantitatively analysed in the 
next subsection. 

C. Time correlation of locations of dissipation 

We now determine the time correlation of the locations of dissipation places. A dissipa- 
tion (a collision) takes place in space r(t) = (x(t),y(t)) at time t. We calculate the power 
spectrum of x{t) and that of y(t) and find that they are of power-law decaying form. Since 
there is no difference between the spectra of x{t) and of y(t), we make average over them 
to obtain P(f). Fig. 4(a) shows the power spectrum for R = 0.1, displaying the form of 
P(f) ~ f~ a with a = 0.49 ± 0.01. The non-zero exponent reflects that clustering occurs 
during the dynamical process: the places where dissipation happens are more likely to be the 
locations for future dissipation. In this figure, the time period corresponding to a frequency 
(/) is t — where 7 is the collision rate in the system. The upper cutoff of frequency is 
fmax = 1024 and the collision rate 7 = 2438. This cutoff of frequecy in Fig. 4(a) corresponds 
to the time interval between two successive collisions (its inverse is the collision rate). Below 
that cutoff, there is no characteristic time scale in the clustering mechanism as indicated 
by the power-law scaling. This spectrum was obtained from the same process as shown in 
Fig. 3(b). We also note that the exponent a depends on the restitution coefficient R. As R 
increases to 1, a decreases to zero. Fig. 4(b) shows the dependence of the exponent a on the 
restitution coefficient. As R increases, the clustering becomes weaker so that the location 
of dissipation places becomes spatially more uniform. This is consistent with the real space 
picture as shown in Fig. 3(b). 
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D. Velocity distribution 



The velocity distribution is found to be Gaussian for all the cases. Fig. 5 shows a typical 
distribution of the square of velocity in a semi-logarithmic plot. The linearity indicates 
that the velocity satisfies the Gaussian distribution like the velocity distribution in usual 
fluids and gases at equilibrium. We will show in Fig. 7(c) that the energy density (e(r) = 
J2iLi ^rnvf5(v — r^)) has no spatial correlation, as expected from the uniform and random 
energy input mechanism. This means that the Gaussian distribution of velocity is valid 
for any part of the system. In such a situation we can now define a granular temperature 
T g =< >, a quantity that is uniform in the system. It should be noted that if the 
velocity does not satisfies the Gaussian distribution the variance of the velocity does not 
have correspondence to the usual temperature but is just internal energy. From the statistical 
physics' point of view, the Gaussian velocity distribution (Maxwellian distribution) comes 
from a very general consideration. In a conservative system, if the momentum enters into 
the Hamiltonian only via the kinetic energy term, the Gaussian distribution of velocity is 



valid irrespective of the potential energy form [^| and the average kinetic energy is ksT 
with ks the Boltzmann constant. 

Non-Gaussian velocity distribution has been found in granular materials under shear 



|23|1 and in fluidized beds of granules 1 . We also find that only by uniformly energy input 



as in the present model without uniaxial force like gravity can the granular materials have 



the Gaussian velocity distribution [^5 



In Fig. 6 we illustrate how the granular temperature T g scales with the energy input rate 
and the inelasticity of the particles. We find 

T g (6V,T,R) = c Q ( x (9) 

where Cq is a constant and ( = ^}J^ ■ Here (5V) 2 /T is the energy input rate. The exponent 
A is found to be 0.65 ± 0.01. We obtain this scaling form from simulations of 75 systems 
with iV = 1024 and p = 0.16 by exploring the three dimensional parameter space (5V, T, R). 
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Here 5V sets the values of 0.001, 0.002, 0.003, 0.005, 0.008, T the values of 0.0001, 0.0005, 
0.001, and R the values of 0.10, 0.45, 0.63, 0.77, 0.89. 



E. Correlations in real space 

We have measured the spatial correlation functions by calculating the following quanti- 
ties: 

Cpp(r) = ^JdR< (p(R, t) - p)(p(R + T ,t)-p)> (10) 
Cv a v P (*) = ^J dR<v a (R,t)vfs(R + r,t) > (11) 

C ee (r) = 1 j dR < (e(R, t) - e) (e(R + r, t) - e) > (12) 

where p(r,t) = S(r — r^t)) is the particle density at position r at time t, 

v a (r) = Y<iLi v ia(t)8( r — r i(t)) is the momentum density in the a-direction and e(r,t) = 
J2iLi \mw^5(r — Yi{t)) is the energy density (here v ia is the a-component of the velocity of 
the i'th particle at time t, a and (3 take the values of x or y and is the position of the i'th 
particle). The < • • • > means time average and the bar means space average. 

From the second rank correlation tensor for the velocity C VaV/3 (r) we define its transverse 
(_L) and logitudinal (||) components 

C\\(r) = r a r/3C VaV0 (r) (13) 

C±(r) =r ±a r ±p C VaV0 (r) (14) 

where r and are unit vectors parallel and perpendicular to the relative displacement r 
respectively. 

For numerical calculations, we first coarse-grain the system into cells of size of about 
5c x 5o". Fig. 7 shows the plots of correlation functions of C pp (r), C\\(r),C±(r) and C ee (r) 
vs. the distance of r. Since the system is isotropic, the correlation functions depend only 
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on the absolute value r = r. From Fig. 7(a), we see clearly the long-range correlation in 
the density fluctuation. For distance less than about 20<r, the correlation is positive while 
for larger distance the correlation is negative. This is in contrast with the situation of 
elastic system of R — 1 where we find C pp {r) is zero for all r > 0. From Fig. 7(b) we 
can see that both parallel and perpendicular components of velocity are correlated in the 
long range. For distance above 15a, the perpendicular velocity correlation is negative. The 



results are different from those for the freely evolving granular medium ]TT[| where there is a 
characteristic length of vortices determined by the negative minimum of the perpendicular 
component of velocity correlation. Here there is no negative minimum in C\\(r) for all r 
less than half of the system size (L ~ 70<r). The lack of a negative minimum indicates no 
characteristic length scale in the system. For the parallel component of velocity correlation, 
one sees that long-range correlation is also built up for the whole system. In comparison, 
we find that for elastic system both parallel and perpendicular components are zero for any 
r > 0. These results make us speculate that our granular system is in a critical state which 
lacks characteristic length in correlations. 

Fig. 7(c) indicates that there are no correlation (too weak to detect, if any) in the energy 
fluctuation, which allows us to name the internal energy as granular temperature in analogy 
to the thermodynamics of usual fluids. 



F. Correlations in reciprocal space 

We have also calculated the correlation in the Fourier space. The structure factors are 
defined as 

S p (k) = i / / < p(r)p(r')e-M r - r ') > drdr' 

(15) 

S v (k) = ±ff< V(r) ■ V(r')e-M r - r ) > drdr' 

where p(r) = J2iLi 5(r — r^) and V(r) = J2iLi Vi<5(r — r^) are the mass density and momentum 
density, and < • • • > means time average. Note that in the above integrals taking r = r' 
will contribute a constant term (irrepective of k), which will be dropped in our calculations 
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in the following, since we are not interested in the self-site correlations. Thus one expects 
that both (S'p(k) and S v (k) for uncorrected systems (e.g. ideal gases) are zero for non-zero 
k. From the simulations with R = 1 (for which the energy-input procedure is not applied), 
we obtain that «SV(k) is indeed zero for non-zero k. However, due to the finite size of the 
particles, small but detectable S p was obtained for non-zero k, which is shown in Fig. 8. 
Let pw\2{r)dV2 denote the probability of finding a particle in dV 2 given that a particle is in 
dVi, S p can be rewritten as 

S p (k) = 1 J dx 2 J d^p\w l2 - l) e *-(«-0 (16) 

Taking the following hard-core exclusion into consideration, 

, r < a 

w 12 = (17) 
I 1 , r > a 

where r = |x 2 — Xi|, it is straightforward to obtain the form of S p (k) 

S p (k)~s + Sl k 2 (for&«i) (18) 

where Sq and si are cr-dependent constants. Fig. 8 indeed shows that 5" p (k) for small values 
of k obeys Eq. flTBp. In the following we present results after subtraction of this finite 
(particle) size effect for the structure factor of mass density. Therefore, by definition, S p (k) 
is zero for non-zero k for elastic systems (R = 1). 

Since the system is isotropic, we obtain the dependence of the structure factors on the 
magnitude of the wave vector k. Let us label the two quantities calculated as described 
above as s p {k) and s v (k). Fig. 9(a) (b) show their dependence on k on log-log plot. The 
linearity indicates that the correlations are of power-law form, indicating no characteristic 
spatial scales in mass density and velocity density fluctuations. The power-law exponent in 
s p (k) ~ fc _/31 is found to keep at a constant value irrespective of R and to our best estimation 
Pi = 1.42 ±0.06 while the exponent in s v (k) ~ k~^ 2 is dependent of the inelastity. Fig. 9(c) 
are the dependence of Pi and P% on the restitution coefficient. One should note that the 
power-law scaling region becomes shorter as R increases. 
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G. Distribution of dissipated energy 



Let us now look at scaling in the energy scale. We focus on all energy dissipation during 
collisions. The energy change after a collision is —5E. Fig. 10(a) shows one time series 
of dissipated energy in the simulation. Fig. 10(b) shows the distributions of 5E counted 
during a long time run for three different restitution coefficients R = 0.6, 0.3 and 0.1. We 
see power-law scaling in the energy scale. The upper cutoff of the scaling (A) is determined 
by the requirement that the energy dissipation rate 7 J A SEP(5E)d(5E) is equal to the 
energy input rate (SV 2 /T) where 7 is the collision rate. We note that the exponent for 
P(8E) ~ 5E~ tl decreases as R increases and the upper cutoff (A) decreases with increasing 
R. For R = 0.1, fi = 0.91 ± 0.02. As R approaches 1, fi decreases to zero. 

H. Anomalous diffusion 

We have traced the trajectories for some randomly marked particles. Fig. 11(a) is a 
typical one. We record the particle positions r(t) = (x(t),y(t)) for a time period of 512 
collisions per particle and calculate the power spectra of x(t) and of y(t) (their average is 
denoted by S(f)). Fig. 11(b) and Fig. 11(c) show S(f) for two different values of R. In 
Fig. 11(b) for R = 0.1, we see that there are two regimes separated by / ~ 100 which 
corresponds to a time period of about 60T where T is the heating period. In the high 
frequency regime, S(f) is proportional to f~ 2 , corresponding to the normal diffusion. Since 
the heating procedure is applied randomly as in Eq. (^), the diffusive behavior is easily 
understood in this short time scale regime. Addtionally, we see from Fig. 11(b) that different 
behavior exists in the lower frequency regime where S(f) is found to be proportional to 
f-(2+p) with p cloge tQ 2 . This gives rise to the mean square displacement 

< (r(t) - r(0)) 2 >~ t l+f} . (19) 

The scaling region of this anomalous diffusion becomes shorter as the restitution coefficient 
increases. This can be seen from Fig. 11(c) which corresponds to R = 0.8. The power-law 
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scaling in the low frequency regime for S(f) is closely related to the power-law scalings 
presented in the above subsections. Since the density fluctuation is correlated in space in a 
self-similar manner as indicated by the power-law scaling in subsection F, the jump of par- 
ticle, depending on the local density fluctuation, may thus satisfy a self-similar distribution. 
This has a link with the Levy-flight random walk 

The anomalous diffusion of Eq. (|i"9l), obtained in two-dimensional simulations, is ap- 
parently similar to the observation by Lewis Fry Richardson in 1926 in the fully developed 
turbulent flow in three dimensions [^f],^]]. Dimensional analysis gives us the correlation of 
the Fourier transformation of velocity 

< v k v. k >~ k~ 5 (20) 

with 5 = (3(3 — l)/(/3 + 1) which is independent of the dimensionality. If (3 = 2, we have 
the Kolmogorov exponent of 5 = 5/3. However, the velocity correlation displayed in Fig. 9 
for R = 0.1 shows that 5 ~ 1.4. At present, we do not know whether this discrepancy is 
intrisic or not. It is possible that an anomalous exponent like the intermittency exponent 
in three-dimensional turbulent flow comes into play. 



IV. DISCUSSION 

The scaling results obtained in the preceding section are robust to the energy-input 
procedure we used. We have verified this by using different parameters for energy-input 
(5V and T) and by using non-periodic energy-input mechanism. We find that the above 
results hold also irrespective of the average density provided that it is in the low density 
regime. For relatively high density, steady state cannot be reached as the particles jam each 
other in such a way that the system is no longer in a "gas" state (i.e., the diffusion constant 
is approaching zero). 

The scaling properties shown in Section III remind us of the concept of self-organized 
criticality (SOC) |2"9"| , |3"0|] . In our model we do not need to fine-tune a control parameter 
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to obtain the scalings in space, time and energy scales. The self-organization appears 
automatically as far as R < 1. However, there are several important differences. In the 
common sense of SOC, avalanches propagate through a medium. The medium is driven 
into a critical state by two opposite processes: external driving and avalanches' disturbance 
(due to threshold instability). In our model, there is no avalanche. Energy is dissipated 
instaneously by collision. Three processes are responsible for the evolution of the system: 
external driving, energy dissipation and self-evolving due to the velocity field. The last 
factor is lacking in the SOC of the sandpile model. 



Our model is also different from turbulence of ordinary fluids plj . In a fully developed 
turbulence, there is a cascade phenomenon which is the origin of the scaling behaviors such 
as the Kolmogorov scaling. The energy is input to the system in a macroscopic scale, and 
then transfered into shorter scales and finally dissipated at a microscopic length scale due to 
viscosity. In our model each particle gains kinetic energy randomly. Therefore, there is no 
macroscopic length (except the system size) asscociated with the energy input. In this sense 
the Reynolds number is negligiblly small in our system. We have obtained an anomalous 
diffusion of a test particle, which is similar to the Richardson law. However, it should be 
noted that our system is in two dimensions and that scaling properties in two dimensions 
are quite different from those in three dimensions in ordinary fluid turbulence. For further 
theoretical studies for the scalings obtained here, one might need to take into account a new 



non-dimensional parameter |32| controlled by dissipation due to inelastic collision. 

In conclusion we have shown by a well-designed model that essential differences exist 
between granular systems and elastic ones. By self-organization our system reaches a critical 
state where there are no characteristic spatial and temporal scales in correlations, and no 
characteristic energy scale in dissipation. 
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FIGURES 

FIG. 1. Energy relaxation starting from a random initial configuration for a system with 
N = 1024, p = 0.16 and R = 0.8. Vertical axis is the kinetic energy per particle (^m < v 2 > with 
m = 2). Best fit using Eq. (7) is also shown with E(oo) = 1.21 x 10" 3 , A = 4.42 x 10 -4 and 
V = 0.38. 

FIG. 2. (a) Real time vs. total number of collisions for system with N = 1024 and p = 0.16. 
The restitution coefficients are R = 0.5, 0.4, 0.3, 0.2 respectively from top to bottom, (b) Collision 
rate (7) vs. degree of dissipation e = 1 — R 2 . The curve is the best fit using Eq. (8) to the 
numerical data with 70 = 505.045, 71 = 130.064 and c = 0.635. 

FIG. 3. (a) Typical configuration of particles. This graph was obtained from simulations for 
system with N = 1024, p = 0.16 and R = 0.2. Arrows represent velocities, (b) Visualization of 
places where dissipation has taken place during a time interval of 20 collisions per particle. When 
a pair of particles collide, we mark a data point (o) at the middle point between the centers of the 
pair. Here N = 1024, R = 0.1, and p = 0.16. 

FIG. 4. (a) Power spectrum P(f) of the time-dependent coordinates of dissipation places 
where collision occurs. The straight line is the best fit to the numerical data, giving a = 0.49 ±0.01. 
(b) Dependence of the exponent a on R. 

FIG. 5. Velocity distribution for system with N = 4096, p = 0.16 and R = 0.6. Vertical axis 
is the number of counts for velocity square in a bin between v 2 and v 2 + Sv 2 with the bin size 
Sv 2 = 1.77 x 10~ 5 . 

FIG. 6. Granular temperature T g vs. ( = v • straight line is the best fit to the 

numerical data, T g = c C A with c = 4.97 x 10 -3 and A = 0.65. 
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FIG. 7. (a) Density-density correlation function C pp ; (b) Parallel (o) and perpendicular 
(+) components of velocity-velocity correlation functions; (c) Energy-energy correlation function. 
Results are obtained for system with N = 1024, p = 0.16 and R = 0.4. Horizontal axis is in unit 
of particle diameter a. 

FIG. 8. Structure factors S p (k) for elastic system (R = 1) of N = 4096 and p = 0.16. Best fit 
using the form of Eq. (18) is also shown with so = 0.55 and s\ = 1.35 x 10~ 6 . 

FIG. 9. (a) Structure factor S p (k); (b) Structure factor S v (k). Power-law decaying exponents 
are ft = 1.39 ± 0.02 and ft = 1.40 ± 0.04 for system with R = 0.1, N = 4096 and p = 0.16. 
Straight lines are best fits to the numerical data in the range of k from 2ir to 100. The length 
scale corresponding to a wave vector k is 2n/k. k=100 corresponds to a length of about 8a. (c) 
Dependence of exponents ft (o) and ft (+) on R. 

FIG. 10. (a) Energy change (-5E) at collision for a process with N = 1024, p = 0.16 and 
R = 0.1. (b) Distribution of dissipated energy P(SE) for three systems with iV = 1024, p = 0.16 
and different R (R = 0.6,0.3,0.1 from right to left). 

FIG. 11. (a) A trajectory for one marked particle taken from simulation with N = 1024, 
p = 0.16 and R = 0.1. (b) Power spectrum S(f) of the time-dependent coordinates of a marked 
particle vs. frequency / for system with iV = 1024, p = 0.16 and R = 0.1. The two lines have 
slopes of —4.17 and —1.97 respectively, (c) Same as (b) but for R = 0.8. The two lines have slopes 
of —4.16 and —1.93 respectively. 
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Fig. 1 -Peng and Ohta 




Fig. 2a -Peng and Ohta 




Fig. 2b -Peng and Ohta 
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Fig. 3a— Peng and Ohta 




Fig. 3b — Peng and Ohta 
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